R Notebook

rm(list = ls())

library(leaflet)
library(maps)
library(rgdal)
library(DT)
library(ggplot2)
library(ggthemes)
library(plotly)
library(magrittr)
library(readxl)
library(plyr)
library(dplyr)
library(tidyr)
library(readr)
library(stringi)
library(RColorBrewer)
library(countrycode)
require(gridExtra)

## IMPORT DATA General ACA information
aca_general_raw <- read_excel("Data/states_aca_general.xlsx", sheet = 2, col_names = FALSE)
# States health outcomes
health_ind_raw <- read_excel("Data/states_health_ind.xlsx", sheet = 4, col_names = FALSE)
# Insurance coverage data Long format
insurance_long <- read_excel("Data/insurance-clean.xlsx", sheet = 1, col_names = TRUE)
# Wide format
insurance <- read_excel("Data/insurance-clean.xlsx", sheet = 2, col_names = TRUE)
# Demographic information
demographics1 <- read.csv("Data/Population_Age_Income.csv", header = TRUE)
# Provided with Assignment 4 (want the yearly population data)
demographics2 <- read_excel("Data/PopulationEstimates.xls")
# Election Results
election <- read_excel("Data/US_Presidential_Results_by_State_1828-2016.xlsx", 
    sheet = 2, col_names = FALSE)

## CLEAN DATA List of states abbreviations
state_abbreviations <- c("AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", 
    "GA", "HI", "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD", "MA", 
    "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ", "NM", "NY", "NC", 
    "ND", "OH", "OK", "OR", "PA", "RI", "SC", "SD", "TN", "TX", "UT", "VT", 
    "VA", "WA", "WV", "WI", "WY")
## CLEAN INSURANCE RATES DATA (WIDE) Drop unnecessary columns
insurance <- insurance[1:50, -1]
# Rename variables
colnames(insurance) <- c("state_abb", "state", "uninsured_num_2008", "uninsured_pct_2008", 
    "uninsured_num_2009", "uninsured_pct_2009", "uninsured_num_2010", "uninsured_pct_2010", 
    "uninsured_num_2011", "uninsured_pct_2011", "uninsured_num_2012", "uninsured_pct_2012", 
    "uninsured_num_2013", "uninsured_pct_2013", "uninsured_num_2014", "uninsured_pct_2014", 
    "uninsured_num_2015", "uninsured_pct_2015")
# Create variables for insurance rates
insurance$insured_pct_2008 <- 100 - insurance$uninsured_pct_2008
insurance$insured_pct_2009 <- 100 - insurance$uninsured_pct_2009
insurance$insured_pct_2010 <- 100 - insurance$uninsured_pct_2010
insurance$insured_pct_2011 <- 100 - insurance$uninsured_pct_2011
insurance$insured_pct_2012 <- 100 - insurance$uninsured_pct_2012
insurance$insured_pct_2013 <- 100 - insurance$uninsured_pct_2013
insurance$insured_pct_2014 <- 100 - insurance$uninsured_pct_2014
insurance$insured_pct_2015 <- 100 - insurance$uninsured_pct_2015
# Reorder variables
insurance <- insurance[, c("state_abb", "state", "uninsured_num_2008", "uninsured_pct_2008", 
    "insured_pct_2008", "uninsured_num_2009", "uninsured_pct_2009", "insured_pct_2009", 
    "uninsured_num_2010", "uninsured_pct_2010", "insured_pct_2010", "uninsured_num_2011", 
    "uninsured_pct_2011", "insured_pct_2011", "uninsured_num_2012", "uninsured_pct_2012", 
    "insured_pct_2012", "uninsured_num_2013", "uninsured_pct_2013", "insured_pct_2013", 
    "uninsured_num_2014", "uninsured_pct_2014", "insured_pct_2014", "uninsured_num_2015", 
    "uninsured_pct_2015", "insured_pct_2015")]
## CLEAN INSURANCE RATES DATA (LONG) Rename variables
colnames(insurance_long) <- c("year", "state", "population", "uninsured_num", 
    "uninsured_pct", "insured_num", "insured_pct")
# Add state abbreviations First create a look-up table (based on wide data)
state_abb_lookup <- insurance[, c("state", "state_abb")]
# Then do the matching
insurance_long$state_abb <- state_abb_lookup$state_abb[match(insurance_long$state, 
    state_abb_lookup$state)]
# Reorder variable
insurance_long <- insurance_long[, c(1:2, ncol(insurance_long), 3:(ncol(insurance_long) - 
    1))]
# Order observations (by state, then by year)
insurance_long <- insurance_long[order(insurance_long$state, insurance_long$year), 
    ]
## CLEAN ACA GENERAL DATA Keep only relevant rows and columns
aca_general <- aca_general_raw[6:57, 1:73]
# Label variables for ACA general data
names(aca_general) <- c("state", "unins_all_pct_10", "unins_all_pct_15", "unins_all_decr_pct", 
    "unins_all_decr", "cov_emp", "cov_parents_plan", "lifetime_lim_preACA_tot", 
    "lifetime_lim_preACA_child", "lifetime_lim_preACA_adultM", "lifetime_lim_preACA_adultF", 
    "cov_private_tot", "cov_private_child", "cov_private_adultM", "cov_private_adultF", 
    "premium_emp_avg_growth_pct_00to10", "premium_emp_avg_growth_pct_10to15", 
    "premium_emp_savings_15", "premium_emp_savings_16", "MLR_rebate_beneficiaries_12", 
    "MLR_rebate_amt_12", "MLR_rebate_beneficiaries_13", "MLR_rebate_amt_13", 
    "MLR_rebate_beneficiaries_14", "MLR_rebate_amt_14", "MLR_rebate_beneficiaries_15", 
    "MLR_rebate_amt_15", "MLR_rebate_amt_12to15", "medicaid_enroll_13", "medicaid_enroll_16", 
    "medicaid_enroll_child_13", "medicaid_enroll_incr_13to16", "medicaid_full_duals", 
    "medicaid_partial_duals", "medicaid_full_or_partial_duals", "state_has_expanded", 
    "insurance_incr_medicaid", "cholest_scr_incr_medicaid", "mammogram_incr_medicaid", 
    "papsmear_incr_medicaid", "clinic_care_incr_medicaid", "all_care_incr_medicaid", 
    "phycisian_visit_ann_incr_medicaid", "depression_decr_medicaid", "good_health_incr_medicaid", 
    "deaths_ann_decr_medicaid", "catastrophic_oop_ann_decr_medicaid", "indebted_ppl_decr_medicaid", 
    "fed_spending_net_incr_inMil", "uncompensated_care_decr_inMil", "mental_substance_elig_medicaid_share", 
    "mental_substance_elig_medicaid", "preexisting_condition_09", "cov_mkt_plan_16", 
    "cov_mkt_16", "receive_tax_credit", "avg_tax_credit", "receive_cost_sharing", 
    "avg_num_mkt_plans_avail_17", "cov_mkt_under75D_pct", "cov_mkt_under100D_pct", 
    "cov_offmkt_elig_tax_credit", "rate_review_funds_to_state", "HIECP_grant_award_to_state", 
    "medicare_enroll_16", "medicare_benef_donuthole", "medicare_benef_donuthole_savings_tot", 
    "medicare_benef_donuthole_savings_avg", "medicare_partB_free_prev_services", 
    "medicare_partB_free_prev_services_share", "medicare_incr_readmit_rate", 
    "medicare_avoided_readmit", "accountable_care_org_num")
# Trim white space on state variable
aca_general$state <- stri_trim_both(aca_general$state)
# Change format of variables Apart from T/F variable for 'state has
# expanded'
aca_general$state_has_expanded[aca_general$state_has_expanded == "yes"] <- TRUE
aca_general$state_has_expanded[aca_general$state_has_expanded == "no"] <- FALSE
# Convert others to numeric form
for (j in c(2:35, 37:73)) {
    aca_general[[j]] <- as.numeric(aca_general[[j]])
}
# Save USA data separately
aca_general_USA <- aca_general[aca_general$state == "United States", ]
# Drop USA and DC data
aca_general <- aca_general[!(aca_general$state %in% c("United States", "District of Columbia")), 
    ]
# Add state abbreviations
aca_general$state_abb <- state_abbreviations
aca_general <- aca_general[, c(1, ncol(aca_general), 2:(ncol(aca_general) - 
    1))]
# Drop raw data
rm(aca_general_raw)
## CLEAN HEALTH INDICATOR DATA Keep only relevant rows and columns
health_ind <- health_ind_raw[4:6703, c(1, 3:6)]
# Label variables
names(health_ind) <- c("state", "measure", "year", "rate", "rank")
# Simplify date (years) - only keep last year if range was given
health_ind$year[health_ind$year == "07/2009 - 06/2012"] <- "2012"
health_ind$year[health_ind$year == "07/2010 - 06/2013"] <- "2013"
health_ind$year[health_ind$year == "07/2012 - 06/2015"] <- "2015"
health_ind$year[health_ind$year == "10/2013-9/2014"] <- "2014"
health_ind$year[health_ind$year == "2008-09"] <- "2009"
health_ind$year[health_ind$year == "2010-11"] <- "2011"
health_ind$year[health_ind$year == "2011-12"] <- "2012"
health_ind$year[health_ind$year == "2011/12"] <- "2012"
health_ind$year[health_ind$year == "2012-13"] <- "2013"
health_ind$year[health_ind$year == "2013-14"] <- "2014"
health_ind$year[health_ind$year == "2013(Q2-Q4)"] <- "2013"
health_ind$year[health_ind$year == "2014-15"] <- "2015"
health_ind$year[health_ind$year == "2015(Q2-Q4)"] <- "2015"
# Export to excel to manually label variable names
unique_measures <- unique(as.factor(health_ind$measure))
# Copied in from excel sheet
unique_measures_relabelled <- c("a.summary_access", "a.unins_adult", "a.unins_child", 
    "a.no_care_bc_cost_adult", "a.high_OOP_relative_under65", "a.at_risk_no_routine_doc_adult", 
    "a.no_dental_adult", "q.summary_prev_treat", "q.with_usual_care_adult", 
    "q.with_cancer_screening_adult", "q.with_vaccines_adult", "q.with_medical_home_child", 
    "q.with_prev_medical_dental_child", "q.with_mental_healthcare_child", "q.with_vaccines_infant", 
    "q.drug_should_avoid_medicare", "q.drug_should_avoid_3conditions_medicare", 
    "q.good_health_provider_medicare", "q.mortality_4conditions_medicare", "q.CLABSI_infection_ratio", 
    "q.info_recovery_hospitalized", "q.good_hospital_staff_hospitalized", "q.improve_mobility_homehealth", 
    "q.improved_wounds_homehealth", "q.sores_NHres", "q.antipsychotic_med_NHres", 
    "u.summary_avoidable_hosp_cost", "u.hosp_asthma_child", "u.hosp_ambulatory_65to74yrs", 
    "u.hosp_ambulatory_above75yrs", "u.30day_hosp_readmit_medicare", "u.30day_hosp_readmit_NHres", 
    "u.hosp_6mos_NHres", "u.hosp_medicare_homehealth", "u.avoidable_ER_medicare", 
    "u.tot_reimb_employer_ins", "u.tot_reimb_medicare", "h.summary_healthy_lives", 
    "h.deaths_amenable", "h.yrs_lost_potential_life_before75", "h.deaths_breast_cancer_F", 
    "h.deaths_colorectal_cancer", "h.deaths_suicide", "h.deaths_infant_mortality", 
    "h.poor_health_adult", "h.smoke_adult", "h.obese_adult", "h.obese_child", 
    "h.poor_dental_adult", "u.premium_emp_private", "u.premium_emp_private_unadj", 
    "u.reimb_medicare_unadj", "u.deaths_amenable_black", "u.deaths_amenable_white", 
    NA, "q.with_prev_screening_above50yrs", NA)
# Replace variable names with abbreviated version
for (i in 1:57) {
    health_ind$measure <- replace(health_ind$measure, health_ind$measure == 
        unique_measures[i], unique_measures_relabelled[i])
}
# Remove the NA's
health_ind <- health_ind[!is.na(health_ind$measure), ]
# Remove duplicates By state, measure.year
health_ind <- health_ind[!duplicated(health_ind[, 1:3]), ]
# Interact measure and year
health_ind$measure_year <- interaction(health_ind$measure, health_ind$year)
# Convert to character variable
health_ind$measure_year <- as.character(health_ind$measure_year)
# Pre-reshape: Store existing data as 'long'
health_ind_long <- health_ind
# Reshape wide
health_ind <- as.data.frame(cbind(health_ind$state, health_ind$measure_year, 
    health_ind$rate))
# Drop rank data for now
names(health_ind) <- c("state", "measure_year", "rate")
health_ind <- reshape(health_ind, idvar = "state", timevar = "measure_year", 
    direction = "wide")
# Drop 'rate' from the name
names(health_ind)[2:ncol(health_ind)] <- substr(names(health_ind), 6, nchar(names(health_ind)))[2:ncol(health_ind)]
# Drop variables who have all entries as NA Due to 'unbalanced panel' in
# reshaping First create a vector to record if variable is all NA
variable_all_NA <- vector(mode = "logical", length = ncol(health_ind))
for (j in 1:ncol(health_ind)) {
    if (sum(is.na(health_ind[[j]])) == 52) {
        variable_all_NA[j] <- TRUE
    }
}
# Only keep the variables that are NOT all NA
health_ind <- health_ind[, !variable_all_NA]
# Change format of variables To character variable
health_ind$state <- as.character(health_ind$state)
# To numeric variable
for (j in c(2:164)) {
    health_ind[[j]] <- as.character(health_ind[[j]])
    health_ind[[j]] <- as.numeric(health_ind[[j]])
}
# Save USA data separately
health_ind_USA <- health_ind[health_ind$state == "United States", ]
# Drop USA and DC data
health_ind <- health_ind[!(health_ind$state %in% c("United States", "District of Columbia")), 
    ]
# Label state abbreviations
health_ind$state_abb <- state_abbreviations
health_ind <- health_ind[, c(1, ncol(health_ind), 2:(ncol(health_ind) - 1))]
# Drop raw data
rm(health_ind_raw)
## CLEAN DEMOGRAPHIC (1) DATA Drop unnecessary row and column
demographics1 <- demographics1[-1, -1]
# Rename variables
colnames(demographics1) <- c("state", "income", "population", "ppl_age0to4", 
    "ppl_age5to9", "ppl_age10to14", "ppl_age15to19", "ppl_age20to24", "ppl_age25to29", 
    "ppl_age30to34", "ppl_age35to39", "ppl_age40to44", "ppl_age45to49", "ppl_age50to54", 
    "ppl_age55to59", "ppl_age60to64", "ppl_age65to69", "ppl_age70to74", "ppl_age75to79", 
    "ppl_age80to84", "ppl_age85plus")
# Drop observations DC and PR
demographics1 <- demographics1[!demographics1$state == "District of Columbia" & 
    !demographics1$state == "Puerto Rico", ]
# Label state abbreviations
demographics1$state_abb <- state_abbreviations
demographics1 <- demographics1[, c(1, ncol(demographics1), 2:(ncol(demographics1) - 
    1))]
# Change format of variables To character variable
demographics1$state <- as.character(demographics1$state)
# To numeric variable
for (j in c(3:ncol(demographics1))) {
    demographics1[[j]] <- as.character(demographics1[[j]])
    demographics1[[j]] <- as.numeric(demographics1[[j]])
}
# BEA regions
demographics1$BEA_region <- ""
demographics1$BEA_region[demographics1$state_abb %in% c("ME", "NH", "VT", "MA", 
    "CT", "RI")] <- "New England"
demographics1$BEA_region[demographics1$state_abb %in% c("NY", "NJ", "PA", "MD", 
    "DE")] <- "Mideast"
demographics1$BEA_region[demographics1$state_abb %in% c("WI", "IL", "MI", "IN", 
    "OH")] <- "Great Lakes"
demographics1$BEA_region[demographics1$state_abb %in% c("WV", "KY", "VA", "TN", 
    "NC", "SC", "AR", "LA", "MS", "AL", "GA", "FL")] <- "Southeast"
demographics1$BEA_region[demographics1$state_abb %in% c("ND", "SD", "NE", "KS", 
    "MN", "IA", "MO")] <- "Plains"
demographics1$BEA_region[demographics1$state_abb %in% c("MT", "ID", "WY", "UT", 
    "CO")] <- "Rocky Mountains"
demographics1$BEA_region[demographics1$state_abb %in% c("AZ", "NM", "TX", "OK")] <- "Southwest"
demographics1$BEA_region[demographics1$state_abb %in% c("WA", "OR", "CA", "NV", 
    "AK", "HI")] <- "Far West"
# Tag income level (category)
demographics1$income_level <- ""
demographics1$income_level[demographics1$income >= 60000] <- "High"
demographics1$income_level[demographics1$income >= 50000 & demographics1$income < 
    60000] <- "Upper Middle"
demographics1$income_level[demographics1$income >= 45000 & demographics1$income < 
    50000] <- "Lower Middle"
demographics1$income_level[demographics1$income < 45000] <- "Low"
# Convert to factor variable
demographics1$income_level <- factor(demographics1$income_level, levels = c("Low", 
    "Lower Middle", "Upper Middle", "High"))
## CLEAN DEMOGRAPHIC (2) DATA
colnames(demographics2) <- demographics2[2, ]
# Drop first few rows (non-observations)
demographics2 <- demographics2[-(1:2), ]
# Drop county-level data, keep state-level data only
demographics2 <- demographics2[demographics2$Area_Name == "United States" | 
    demographics2$State != lag(demographics2$State), ]
# Keep only columns of population from 2010 to 2015
demographics2 <- demographics2[, c("State", "Area_Name", "POP_ESTIMATE_2010", 
    "POP_ESTIMATE_2011", "POP_ESTIMATE_2012", "POP_ESTIMATE_2013", "POP_ESTIMATE_2014", 
    "POP_ESTIMATE_2015")]
# Change format of variables To character variable
for (j in c(3:8)) {
    demographics2[[j]] <- as.numeric(demographics2[[j]])
}
# Rename variables
colnames(demographics2) <- c("state_abb", "state", "population_2010", "population_2011", 
    "population_2012", "population_2013", "population_2014", "population_2015")
## CLEAN ELECTION DATA Keep only 2012 and 2008 data, and data for the 50
## states
election <- election[4:53, c("X0", "X5", "X6", "X9", "X10")]
# Rename variables
colnames(election) <- c("state", "Dem_pct_2012", "Rep_pct_2012", "Dem_pct_2008", 
    "Rep_pct_2008")
# Add state abbreviations Match with look-up table
election$state_abb <- state_abb_lookup$state_abb[match(election$state, state_abb_lookup$state)]
# Change format of variables To numeric variable
for (j in c(4:6)) {
    election[[j]] <- as.numeric(election[[j]])
}
# Calculate winning party for each state
election$party_2012 <- ifelse(election$Dem_pct_2012 > election$Rep_pct_2012, 
    "Democratic", "Republican")
election$party_2008 <- ifelse(election$Dem_pct_2008 > election$Rep_pct_2008, 
    "Democratic", "Republican")
# Reorder variables
election <- election[, c("state", "state_abb", "Dem_pct_2012", "Rep_pct_2012", 
    "party_2012", "Dem_pct_2008", "Rep_pct_2008", "party_2008")]

## EXPORT DATASETS FOR SHINY USE
#write_csv(insurance, "/Users/gracekongyx/Documents/*6_Data Science/QMSS G4063 Data Visualization/Final Project/Data/shiny_insurance.csv")

(1) Health Insurance Coverage

## SOME DATA PREPARATION Look-up information into insurance and
## insurance_long, on: Party (2012)
insurance_long$party_2012 <- election$party_2012[match(insurance_long$state, 
    election$state)]
insurance$party_2012 <- election$party_2012[match(insurance$state, election$state)]
# BEA region
insurance_long$BEA_region <- demographics1$BEA_region[match(insurance_long$state, 
    demographics1$state)]
insurance$BEA_region <- demographics1$BEA_region[match(insurance$state, demographics1$state)]
# Income level
insurance_long$income_level <- demographics1$income_level[match(insurance_long$state, 
    demographics1$state)]
insurance$income_level <- demographics1$income_level[match(insurance$state, 
    demographics1$state)]
# Also find nation-wide trend for insurance
insurance_long_US <- data.frame(year = c(2008, 2009, 2010, 2011, 2012, 2013, 
    2014, 2015), uninsured_pct = c(16.8, 17.5, 18.2, 17.2, 16.9, 16.6, 13.3, 
    10.5))

The following figure displays the trend in the uninsured rate between 2008 and 2015. Results are shown for the USA average (in black), as well as for every individual state (in color). The states are grouped and color-coded by region to facilitate identification, where one can toggle the display to show a particular region to locate a particular state. We see a clear downward trend in the uninsured rate once most of the Affordable Care Act (ACA) provisions were implemented in 2014.

## PLOT 1
plot1 <- (ggplot(insurance_long, aes(year, uninsured_pct)) + geom_line(aes(color = BEA_region, 
    group = state, label = uninsured_num)) + scale_x_continuous(breaks = c(2008, 
    2009, 2010, 2011, 2012, 2013, 2014, 2015)) + scale_y_continuous(breaks = c(5, 
    10, 15, 20, 25, 30)) + labs(x = "Year", y = "% Uninsured", color = "Region") + 
    geom_line(data = insurance_long_US, lwd = 0.8) + geom_vline(xintercept = 2010, 
    linetype = 2, color = "gray30") + geom_vline(xintercept = 2014, linetype = 2, 
    color = "gray30") + annotate("text", x = 2009.4, y = 25, label = "Obama \nsigns ACA", 
    color = "black", fontface = 2, size = 3) + annotate("text", x = 2013.4, 
    y = 24.5, label = "Most ACA \nprovisions \ntake effect", color = "black", 
    fontface = 2, size = 3) + annotate("text", x = 2014.2, y = 11, label = "USA \nAverage", 
    color = "black", fontface = 2, size = 3) + ggtitle("Change in Percentage Uninsured by State, \nfrom 2008 to 2015") + 
    theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0, 
    size = 12)))
ggplotly(plot1, dynamicTicks = FALSE, tooltip = c("state", "year", "uninsured_pct", 
    "uninsured_num"))

In the following plot, I explore the difference between Democratic and Republican states, based on the 2012 US Presidential Election results. It is clear that the democratic states largely have a lower uninsured rate than the Republican states, though both groups experienced declines in the uninsured rates in 2014.

## PLOT 2
plot2 <- (ggplot(insurance_long, aes(year, uninsured_pct)) + geom_line(aes(color = party_2012, 
    group = state, label = uninsured_num)) + scale_x_continuous(breaks = c(2008, 
    2009, 2010, 2011, 2012, 2013, 2014, 2015)) + scale_y_continuous(breaks = c(5, 
    10, 15, 20, 25, 30)) + scale_color_manual(values = c("steel blue", "firebrick1")) + 
    labs(x = "Year", y = "% Uninsured", color = "Party \n(2012 Winner)") + geom_line(data = insurance_long_US, 
    lwd = 1) + geom_vline(xintercept = 2010, linetype = 2, color = "gray30") + 
    geom_vline(xintercept = 2014, linetype = 2, color = "gray30") + annotate("text", 
    x = 2009.4, y = 25, label = "Obama \nsigns ACA", color = "black", fontface = 2, 
    size = 3) + annotate("text", x = 2013.4, y = 24.5, label = "Most ACA \nprovisions \ntake effect", 
    color = "black", fontface = 2, size = 3) + annotate("text", x = 2014.2, 
    y = 11, label = "USA \nAverage", color = "black", fontface = 2, size = 3) + 
    ggtitle("Change in Percentage Uninsured by Party, \nfrom 2008 to 2015") + 
    theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0, 
    size = 12)))
ggplotly(plot2, dynamicTicks = FALSE, tooltip = c("state", "year", "uninsured_pct", 
    "uninsured_num"))

## PLOT 3 Subset relevant variables
insurance_2015 <- insurance[, c("state", "uninsured_pct_2015", "insured_pct_2015", 
    "party_2012", "BEA_region", "income_level")]
# Boxplots by party
plot3a <- (ggplot(insurance_2015, aes(x = reorder(party_2012, insured_pct_2015), 
    y = insured_pct_2015)) + geom_boxplot(aes(fill = party_2012, label = insured_pct_2015), 
    color = "gray30") + coord_flip() + scale_fill_manual(values = c("steel blue", 
    "firebrick")) + ggtitle("Democratic vs. Republican States") + theme_minimal() + 
    theme(plot.title = element_text(face = "bold", hjust = 0.37, size = 12), 
        legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank()))
# Boxplots by BEA region
plot3b <- (ggplot(insurance_2015, aes(x = reorder(BEA_region, insured_pct_2015), 
    y = insured_pct_2015)) + geom_boxplot(aes(fill = BEA_region, label = insured_pct_2015), 
    color = "gray30") + coord_flip() + scale_fill_manual(values = brewer.pal(n = 8, 
    name = "Dark2")) + ggtitle("By Bureau of Economic Analysis (BEA) Region") + 
    theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.24, 
    size = 12), legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank()))
# Boxplots by income level
plot3c <- (ggplot(insurance_2015, aes(x = income_level, y = insured_pct_2015)) + 
    geom_boxplot(aes(fill = income_level, label = insured_pct_2015), color = "gray30") + 
    coord_flip() + scale_fill_manual(values = brewer.pal(n = 4, name = "Dark2")) + 
    labs(x = "", y = "% Insured (2015)") + ggtitle("By State's Income Level") + 
    theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.39, 
    size = 12), legend.position = "none", axis.title.y = element_blank()))
# Combine plots vertically
grid.arrange(plot3a, plot3b, plot3c, nrow = 3, top = "2015 Insurance Rates at a Glance")

## SOME DATA RESHAPING Reshape aca_general (long), keeping the interesting
## quantities Import uninsured number (2015)
aca_general$unins_all_15 <- insurance$uninsured_num_2015[match(aca_general$state, 
    insurance$state)]
# Keep only select variables
aca_general_long <- aca_general[, c("state", "state_abb", "cov_emp", "cov_mkt_plan_16", 
    "medicaid_enroll_16", "medicare_enroll_16", "unins_all_15")]
# The actual reshaping
aca_general_long <- reshape(aca_general_long, varying = c("cov_emp", "cov_mkt_plan_16", 
    "medicaid_enroll_16", "medicare_enroll_16", "unins_all_15"), v.names = "Number", 
    timevar = "Source", times = c("Employer", "Marketplace", "Medicaid / CHIP", 
        "Medicare", "Uninsured"), new.row.names = 1:1000, direction = "long")
# Look-up information into aca_general_long, on: Party (2012)
aca_general_long$party_2012 <- election$party_2012[match(aca_general_long$state, 
    election$state)]
# BEA region
aca_general_long$BEA_region <- demographics1$BEA_region[match(aca_general_long$state, 
    demographics1$state)]
# Income level
aca_general_long$income_level <- demographics1$income_level[match(aca_general_long$state, 
    demographics1$state)]
# Population
aca_general_long$population_2015 <- demographics2$population_2015[match(aca_general_long$state, 
    demographics2$state)]
# Calculate percentage from number and population
aca_general_long$percentage <- aca_general_long$Number/aca_general_long$population_2015 * 
    100

Next, I investigate the main health insurance and coverage sources, namely employer coverage, the individual marketplace, Medicaid / Childrens’ Health Insurance Program (CHIP), and Medicare for the aged. I broke down the results by party and income level, taking a simple (unweighted) average across the states in the categories. We observe that employer coverage rates increase as income level increases, especially for Democratic states. Democratic and low income states have more people on Medicaid / CHIP than Republican and low income states.

## PLOT 4
aca_general_long_by_party_income <- ddply(aca_general_long, c("party_2012", 
    "income_level", "Source"), summarise, avg_pct = mean(percentage))
plot4 <- (ggplot(aca_general_long_by_party_income, aes(income_level, avg_pct, 
    fill = Source)) + geom_col(position = "dodge") + facet_grid(party_2012 ~ 
    .) + labs(x = "Income Level", y = "Average % in State", color = "Source of \nInsurance") + 
    scale_fill_manual(values = brewer.pal(n = 6, name = "Dark2")[c(1:3, 6, 4)]) + 
    ggtitle("Sources of Health Insurance, \nby Party and Income Level") + theme_minimal() + 
    theme(plot.title = element_text(face = "bold", hjust = 0, size = 12)))
ggplotly(plot4, tooltip = c("income_level", "Source", "avg_pct"))

## PREPARATION FOR FURTHER PLOTS (INCLUDING DATA TABLE) Want to combine
## information on insurance coverage with key information on health
## outcomes/access and affordability
key_indicators <- state_abb_lookup
## Merge in information to the data table on the following Demographic
## Information Income level
key_indicators$income_level <- demographics1$income_level[match(key_indicators$state_abb, 
    demographics1$state_abb)]
# Insurance Coverage Percentage insured (2015)
key_indicators$insured_pct_2015 <- insurance$insured_pct_2015[match(key_indicators$state_abb, 
    insurance$state_abb)]
# Percentage insured (2014)
key_indicators$insured_pct_2014 <- insurance$insured_pct_2014[match(key_indicators$state_abb, 
    insurance$state_abb)]
# Indicator: Whether state has expanded medicare
key_indicators$state_has_expanded <- aca_general$state_has_expanded[match(key_indicators$state_abb, 
    aca_general$state_abb)]
# Health outcomes and access Mortality amenable to health care, deaths per
# 10000 population
key_indicators$deaths_amenable_2014 <- health_ind$h.deaths_amenable.2014[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Years of Potential Life Lost before 75
key_indicators$years_life_lost <- health_ind$h.yrs_lost_potential_life_before75.2014[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Adults with a usual source of care (%)
key_indicators$usual_care_2015 <- health_ind$q.with_usual_care_adult.2015[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Adults with age/gender-appropriate cancer screenings (%)
key_indicators$with_cancer_screening <- health_ind$q.with_cancer_screening_adult.2014[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Adults with age-appropriate vaccines (%)
key_indicators$with_vaccines <- health_ind$q.with_vaccines_adult.2015[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Affordability and cost-efficiency Individuals under age 65 with high OOP
# medical costs relative to annual household income
key_indicators$high_OOP_relative <- health_ind$a.high_OOP_relative_under65.2015[match(key_indicators$state_abb, 
    health_ind$state_abb)]
# Average annual growth in family premiums for employer coverage (between
# 2010 and 2015)
key_indicators$premium_ann_growth_10_15 <- aca_general$premium_emp_avg_growth_pct_10to15[match(key_indicators$state_abb, 
    aca_general$state_abb)]
# Marketplace consumers who could select a plan for less than $100
key_indicators$IM_plan_under_100 <- aca_general$cov_mkt_under100D_pct[match(key_indicators$state_abb, 
    aca_general$state_abb)]
# Net increase in federal spending (millions)
key_indicators$incr_fed_spending_mil <- aca_general$fed_spending_net_incr_inMil[match(key_indicators$state_abb, 
    aca_general$state_abb)]

(2) Health Outcomes and Access

## PLOT 5
plot5 <- (ggplot(key_indicators, aes(insured_pct_2014, deaths_amenable_2014, 
    color = income_level, shape = state_has_expanded)) + geom_point(size = 2) + 
    geom_smooth(aes(group = income_level), method = "lm", lwd = 1, se = FALSE) + 
    labs(x = "% Insured", y = "Deaths Preventable by Healthcare Intervention \n(per 100,000 People)", 
        color = "Income \nLevel", shape = "State has \nExpanded \nMedicaid") + 
    scale_color_manual(values = brewer.pal(n = 4, name = "Dark2")) + ggtitle("Relationship Between Avoidable Mortality and \nHealth Insurance Coverage Among States") + 
    theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5, 
    size = 12)))
plot5

## DATA PROCESSING Obtain shape files of US states
map_states = map("state", fill = TRUE, plot = FALSE)
# Standardize name format with rest of assignment Function to conver to
# proper case
properCase <- function(x) {
    s <- strsplit(x, " ")[[1]]
    paste(toupper(substring(s, 1, 1)), substring(s, 2), sep = "", collapse = " ")
}
# Convert to proper case
map_states$state <- sapply(map_states$names, properCase)
# Rename the main component of states with multiple parts in the map
map_states$state[map_states$state == "Massachusetts:main"] <- "Massachusetts"
map_states$state[map_states$state == "Michigan:south"] <- "Michigan"
map_states$state[map_states$state == "New York:main"] <- "New York"
map_states$state[map_states$state == "North Carolina:main"] <- "North Carolina"
map_states$state[map_states$state == "Virginia:main"] <- "Virginia"
map_states$state[map_states$state == "Washington:main"] <- "Washington"
# Match with state abbreviations
map_states$state_abb <- state_abb_lookup$state_abb[match(map_states$state, state_abb_lookup$state)]
# Then strip off all the extra location information for states with multiple
# parts (leaving just the state name)
for (i in 1:length(map_states$state)) {
    map_states$state[i] <- unlist(strsplit(map_states$state[i], ":"))[1]
}
# Get coordinates of state centers
states_centers <- state.center
state.center <- cbind(states_centers, state_abb_lookup)
# Health outcomes: Import key indicators data, matched by state Remember to
# later match information based on state, not state abbreviation (due to how
# we labelled above)
map_states$deaths_amenable_2014 <- key_indicators$deaths_amenable_2014[match(map_states$state, 
    key_indicators$state)]
map_states$usual_care_2015 <- key_indicators$usual_care_2015[match(map_states$state, 
    key_indicators$state)]
map_states$with_cancer_screening <- key_indicators$with_cancer_screening[match(map_states$state, 
    key_indicators$state)]
map_states$with_vaccines <- key_indicators$with_vaccines[match(map_states$state, 
    key_indicators$state)]
# Cost outcomes: Import key indicators data, matched by state Remember to
# later match information based on state, not state abbreviation (due to how
# we labelled above)
map_states$high_OOP_relative <- key_indicators$high_OOP_relative[match(map_states$state, 
    key_indicators$state)]
map_states$premium_ann_growth_10_15 <- key_indicators$premium_ann_growth_10_15[match(map_states$state, 
    key_indicators$state)]
map_states$IM_plan_under_100 <- key_indicators$IM_plan_under_100[match(map_states$state, 
    key_indicators$state)]
map_states$incr_fed_spending_mil <- key_indicators$incr_fed_spending_mil[match(map_states$state, 
    key_indicators$state)]

## PLOT 6: LEAFLET MAP - HEALTH OUTCOMES
(leaflet(map_states) %>% setView(lat = 39.8282, lng = -96, zoom = 4) %>% addPolygons(color = "#333333", 
    weight = 1, smoothFactor = 0.5, fillOpacity = 0) %>% addPolygons(group = "Preventable Deaths", 
    fillColor = ~colorQuantile("RdYlGn", -deaths_amenable_2014)(-deaths_amenable_2014), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Preventable Mortality</b>", "<br/>", 
        "<b>per 100,000: </b>", map_states$deaths_amenable_2014)) %>% addPolygons(group = "Access to Usual Care", 
    fillColor = ~colorQuantile("RdYlGn", usual_care_2015)(usual_care_2015), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Access to Usual</b>", "<br/>", "<b>Source of Care: </b>", 
        map_states$usual_care_2015, "%")) %>% addPolygons(group = "Cancer Screenings Rate", 
    fillColor = ~colorQuantile("RdYlGn", with_cancer_screening)(with_cancer_screening), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Cancer Screenings</b>", "<br/>", "<b>Rate: </b>", 
        map_states$with_cancer_screening, "%")) %>% addPolygons(group = "Adult Vaccination Rate", 
    fillColor = ~colorQuantile("RdYlGn", with_vaccines)(with_vaccines), smoothFactor = 0.5, 
    stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", map_states$state, 
        "<br/>", "<b>Vaccination</b>", "<br/>", "<b>Rate: </b>", map_states$with_vaccines, 
        "%")) %>% addLabelOnlyMarkers(data = filter(state.center, state_abb != 
    "AK" & state_abb != "HI"), lng = ~x, lat = ~y, label = ~state_abb, labelOptions = labelOptions(noHide = T, 
    direction = "top", textOnly = T)) %>% addLayersControl(baseGroups = c("Preventable Deaths", 
    "Access to Usual Care", "Cancer Screenings Rate", "Adult Vaccination Rate"), 
    options = layersControlOptions(collapsed = FALSE)))

(3) Affordability and Cost Outcomes

## PLOT 7: LEAFLET MAP - AFFORDABILITY / COST OUTCOMES
(leaflet(map_states) %>% setView(lat = 39.8282, lng = -96, zoom = 4) %>% addPolygons(color = "#333333", 
    weight = 1, smoothFactor = 0.5, fillOpacity = 0) %>% addPolygons(group = "Contained OOP Costs", 
    fillColor = ~colorQuantile("RdYlGn", -high_OOP_relative)(-high_OOP_relative), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>High OOP Relative</b>", "<br/>", "<b>to Income: </b>", 
        map_states$high_OOP_relative, "%")) %>% addPolygons(group = "Premium Growth Rate", 
    fillColor = ~colorQuantile("RdYlGn", -premium_ann_growth_10_15)(-premium_ann_growth_10_15), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Annual Premium</b>", "<br/>", "<b>Growth Rate: </b>", 
        map_states$premium_ann_growth_10_15, "%")) %>% addPolygons(group = "Affordable Individual Marketplace Plan", 
    fillColor = ~colorQuantile("RdYlGn", IM_plan_under_100)(IM_plan_under_100), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Marketplace Consumers</b>", "<br/>", 
        "<b>Who Can Select Plan</b>", "<br/>", "<b>Under $100: </b>", map_states$IM_plan_under_100 * 
            100, "%")) %>% addPolygons(group = "Increase in Federal Spending", 
    fillColor = ~colorQuantile("RdYlGn", -incr_fed_spending_mil)(-incr_fed_spending_mil), 
    smoothFactor = 0.5, stroke = FALSE, fillOpacity = 0.6, popup = paste("<b>State: </b>", 
        map_states$state, "<br/>", "<b>Net Increase in Federal</b>", "<br/>", 
        "<b>Spending</b>: $", map_states$incr_fed_spending_mil, "million")) %>% 
    addLabelOnlyMarkers(data = filter(state.center, state_abb != "AK" & state_abb != 
        "HI"), lng = ~x, lat = ~y, label = ~state_abb, labelOptions = labelOptions(noHide = T, 
        direction = "top", textOnly = T)) %>% addLayersControl(baseGroups = c("Contained OOP Costs", 
    "Premium Growth Rate", "Affordable Individual Marketplace Plan", "Increase in Federal Spending"), 
    options = layersControlOptions(collapsed = FALSE)))

(4) Summary of States’ Performance

# Create data table reporting values
data_table_values <- key_indicators[ , c("state", "insured_pct_2015", "deaths_amenable_2014", "usual_care_2015", "high_OOP_relative", "premium_ann_growth_10_15")]
# Create data table reporting rank
data_table_rank <- data_table_values
# Convert each function to rank
  # Rank 1 is always more favorable (regardless of definition of variable)    
data_table_rank$r_insured_pct_2015 <- rank(-data_table_values$insured_pct_2015, na.last = "keep", ties.method = "min")
data_table_rank$r_deaths_amenable_2014 <- rank(data_table_values$deaths_amenable_2014, na.last = "keep", ties.method = "min")
data_table_rank$r_usual_care_2015 <- rank(-data_table_values$usual_care_2015, na.last = "keep", ties.method = "min")
data_table_rank$r_high_OOP_relative <- rank(data_table_values$high_OOP_relative, na.last = "keep", ties.method = "min")
data_table_rank$r_premium_ann_growth_10_15 <- rank(data_table_values$premium_ann_growth_10_15, na.last = "keep", ties.method = "min")
  # Remove value variables
data_table_rank$insured_pct_2015 <- NULL
data_table_rank$deaths_amenable_2014 <- NULL
data_table_rank$usual_care_2015 <- NULL
data_table_rank$high_OOP_relative <- NULL
data_table_rank$premium_ann_growth_10_15 <- NULL

The following data table presents the rankings of the states on the following insurance coverage, health and affordability outcomes: 1) Percentage of population insured 2) Deaths preventable by healthcare intervention (per 100,000 people) 3) Proportion of adults with access to a usual source of care 4) Proportion of people with high out-of-pocket costs (relative to household income) 5) Annual rate of health insurance premium growth from 2010 to 2015 A rank of 1 always indicates the more favorable outcome (i.e. higher insurance and treatment access rate, lower preventable mortality, OOP costs and premium growth)

## PLOT 8: DATA TABLE RANKING THE STATES ON VARIOUS MEASURES
datatable(data_table_rank, rownames = FALSE, colnames = c("State", "Percentage Insured", 
    "Fewer Preventable Deaths", "Usual Care Access", "OOP Contained", "Lower Annual Premium Growth"), 
    caption = "Rank of the 50 States in Insurance Coverage and Health Outcomes")

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

library(tm)
library(tm.plugin.lexisnexis)
library(readxl)
library(gtools)  # for smartbind
library(dplyr)   # for data_frame
library(lubridate)   # for date formatting
library(stringr)  
library(tools)  # Title case
library(quanteda)
library(ggplot2)

# Combine CSV and HTML Files
data <- read_excel("LexisNexis/NYTimes_Metadata.xlsx")
colnames(data) <- tolower(colnames(data))

# Correct data
data$date <- substr(parse_date_time(data$date, c("mdy")),1,10)
data$author <- toTitleCase(tolower(data$byline))
data$byline <- NULL

## Get Text files
source1 <- LexisNexisSource("LexisNexis/The_New_York1.html")
source2 <- LexisNexisSource("LexisNexis/The_New_York2.html")

corpus1 <- Corpus(source1, readerControl = list(language = NA))
corpus2 <- Corpus(source2, readerControl = list(language = NA))

corpus <- c(corpus1, corpus2)

# Convert to quanteda corpus
corpus <- quanteda::corpus(corpus)

## Add Metadata
# Check: match(data$headline, corpus$documents$heading)
corpus$documents$datetimestamp <- substring(corpus$documents$datetimestamp, 1,4)
corpus$documents$date <- corpus$documents$datetimestamp
corpus$documents$description <- corpus$documents$id

# options(width = 200)
# kwic(corpus, "Trump")
# 
save(corpus, file="nytimes.rda")

the Nummber of articles mentions Affordable care by years

dat <- data.frame(date =corpus$documents$datetimestamp, levels=c(1:1000))
ggplot(data=dat, aes(x=date))+geom_bar(stat="count")+ggtitle("Number of articles mentioned ACA by years")+theme_tufte()+ylab("Number of articles") + xlab("year")

load("nytimes.rda")
dfmtotal <- dfm(corpus, remove = stopwords("english"), stem = TRUE, removePunct = TRUE, removeNumbers = TRUE, tolower = TRUE, verbose = TRUE)

Creating a dfm from a corpus ...
   ... lowercasing
   ... tokenizing
   ... found 1,000 documents, 22,386 features
   ... removed 169 features, from 174 supplied (glob) feature types
   ... stemming features (English), trimmed 7834 feature variants
   ... created a 1,000 x 14,383 sparse dfm
   ... complete. 
Elapsed time: 2.81 seconds.

dfmtotal[, 1:5]

Document-feature matrix of: 1,000 documents, 5 features (88.9% sparse).

head(stopwords("english"), 20)

 [1] "i"          "me"         "my"         "myself"     "we"         "our"        "ours"       "ourselves" 
 [9] "you"        "your"       "yours"      "yourself"   "yourselves" "he"         "him"        "his"       
[17] "himself"    "she"        "her"        "hers"      

freq<-topfeatures(dfmtotal, 20)
library(ggplot2)  
wf <- data.frame(word=names(freq), freq=freq)   
p <- ggplot(subset(wf, freq>50), aes(word, freq))    
p <- p + geom_bar(stat="identity")
p <- p + geom_text(aes(label=freq), vjust=-0.2)+ scale_x_discrete(limits= wf$word)   
p <- p + theme(axis.text.x=element_text(angle=45, hjust=1))+ggtitle("Word frequency from 2011-2017")+theme_tufte()+ylab("Frequency of Stem") + xlab("Top 30 Used Stem")
p

library(quanteda)
library(stringr)
load('nytimes.rda')
library(tm)
NYT_source <- VectorSource(corpus$documents[,1])
NYT <- VCorpus(NYT_source)
documents <- corpus$documents
library(qdap)
removeNumPunct <- function(x){gsub("[^[:alpha:][:space:]]*", "", x)}
clean_corpus <- function(corpus){
corpus <- tm_map(corpus, removePunctuation)
corpus <- tm_map(corpus, content_transformer(tolower))
corpus <- tm_map(corpus, content_transformer(replace_symbol))
corpus <- tm_map(corpus, removeWords, c(stopwords("en"), "will", "can", "Paul
Krugman", "David Brooks", "Thomas L. Friedman", "Maureen Dowd", "Nicholas
Kristof"))
# We could add more stop words as above
corpus <- tm_map(corpus, stripWhitespace)
corpus <- tm_map(corpus, removeNumbers)
corpus <- tm_map(corpus, content_transformer(removeNumPunct))
return(corpus)
}
library(SnowballC)
NYT_clean <- clean_corpus(NYT)
NYT_stem <- tm_map(NYT_clean, stemDocument)
meta(NYT_stem, type = "local", tag = "author") <- documents$author
NYT_dtm <- DocumentTermMatrix(NYT_stem)
NYT_tdm <- TermDocumentMatrix(NYT_stem)
library(dplyr)
library(tidytext)
NYT_tdm2 <- tidy(NYT_tdm)
NYT_dtm2 <- tidy(NYT_dtm)
NYT_tidy <- tidy(NYT_stem)
NYT_tdm2 <- merge(NYT_tdm2, NYT_tidy, by.x = "document", by.y = "id", all.x =
TRUE)
NYT_dtm2 <- merge(NYT_dtm2, NYT_tidy, by.x = "document", by.y = "id", all.x =
TRUE)

library(wordcloud)
atxt <- documents$texts[documents$datetimestamp == "2011"]
btxt <- documents$texts[documents$datetimestamp == "2012"]
ctxt <- documents$texts[documents$datetimestamp == "2013"]
dtxt <- documents$texts[documents$datetimestamp == "2014"]
etxt <- documents$texts[documents$datetimestamp == "2015"]
ftxt <- documents$texts[documents$datetimestamp == "2016"]
gtxt <- documents$texts[documents$datetimestamp == "2017"]
clean.text <- function(x)
{
# tolower
x = tolower(x)
# remove rt
x = gsub("rt", "", x)
# remove at
x = gsub("@\\w+", "", x)
# remove punctuation
x = gsub("[[:punct:]]", "", x)
# remove numbers
x = gsub("[[:digit:]]", "", x)
# remove links http
x = gsub("http\\w+", "", x)
# remove tabs
x = gsub("[ |\t]{2,}", "", x)
# remove blank spaces at the beginning
x = gsub("^ ", "", x)
# remove blank spaces at the end
x = gsub(" $", "", x)
return(x)
}
aclean <- clean.text(atxt)
bclean <- clean.text(btxt)
cclean <- clean.text(ctxt)
dclean <- clean.text(dtxt)
eclean <- clean.text(etxt)
fclean <- clean.text(ftxt)
gclean <- clean.text(gtxt)
a <- paste(aclean, collapse=" ")
b <- paste(bclean, collapse=" ")
c <- paste(cclean, collapse=" ")
d <- paste(dclean, collapse=" ")
e <- paste(eclean, collapse=" ")
f <- paste(fclean, collapse=" ")
g <- paste(gclean, collapse=" ")
# put everything in a single vector
all <- c(a,b,c,d,e,f,g)
all <- removeWords(all, c("will", "can", "the", "that","are"))
# create corpus
corpus2 <- Corpus(VectorSource(all))
# create term-document matrix
cloudtdm <- TermDocumentMatrix(corpus2)
# convert as matrix
cloudtdm <- as.matrix(cloudtdm)
# add column names
colnames(cloudtdm) <- c("2011","2012","2013","2014","2015","2016","2017")
comparison.cloud(cloudtdm, random.order=FALSE,colors = brewer.pal(8, "Dark2"),title.size=1.5, max.words = 200)